;;;; 
;;;; least-squares.lisp
;;;;
;;;; Least Squares routines for linear systems defined
;;;; by Ax=b

(in-package :lib.gl)

(defun column-space-projection (a b)
  (m* a (inverse (m* (transpose a) a))
      (transpose a)
      b))

(defun projection-matrix (a)
  (m* a (inverse (m* (transpose a) a))
      (transpose a)))


(defun linear-least-squares (x-values y-values)
  (let ((a (make-matrix (length x-values) 2 1.0d0))
	(b (make-matrix (length y-values) 1)))
    (set-column a 1 (apply #'%gl-vector x-values))
    (set-column b 0 (apply #'%gl-vector y-values))
    (print-matrix (inverse (m* (transpose a) a)))
    (let* ((at (transpose a))
	   (aainv (inverse (m* at a))))
      (m* aainv (m* at b)))))

;;; Tests

(defun test-column-projection-1 ()
  (assert (matrix-equal  (column-space-projection (%matrix '((1 2)
							     (1 3)
							     (0 0)))
						  (transpose (%matrix '((4 5 6)))))
			 (transpose (%matrix '((4 5 0)))))))